A numerical study of spatio-temporal COVID-19 vaccine model via finite-difference operator-splitting and meshless techniques

In this paper, a new spatio-temporal model is formulated to study the spread of coronavirus infection (COVID-19) in a spatially heterogeneous environment with the impact of vaccination. Initially, a detailed qualitative analysis of the spatio-temporal model is presented. The existence, uniqueness, positivity, and boundedness of the model solution are investigated. Local asymptotical stability of the diffusive COVID-19 model at steady state is carried out using well-known criteria. Moreover, a suitable nonlinear Lyapunov functional is constructed for the global asymptotical stability of the spatio-temporal model. Further, the model is solved numerically based on uniform and non-uniform initial conditions. Two different numerical schemes named: finite difference operator-splitting and mesh-free operator-splitting based on multi-quadratic radial basis functions are implemented in the numerical study. The impact of diffusion as well as some pharmaceutical and non-pharmaceutical control measures, i.e., reducing an effective contact causing infection transmission, vaccination rate and vaccine waning rate on the disease dynamics is presented in a spatially heterogeneous environment. Furthermore, the impact of the aforementioned interventions is investigated with and without diffusion on the incidence of disease. The simulation results conclude that the random motion of individuals has a significant impact on the disease dynamics and helps in setting a better control strategy for disease eradication.

Moreover, the susceptible individuals get infected after an interaction at the rate β and βψ ν with infected and exposed individuals respectively. The fractions of infected individuals among total population are IS N and ES N respectively and therefore, force of infection is, The detailed description of the model (1) embedded parameters along with their values are presented in Table 1. The constants α S , α E , α I , α V and α R represents the diffusion rate of respective populations. Moreover, the following no flux boundary conditions are taken into account for the above model: (x, t) = β (I(x, t) + ψE(x, t))S(x, t) N(x, t) . (2) Initial conditions. The following set of initial conditions, provided by (3) and (4) are utilized to simulate the model (1). Based on the conditions listed in 6 , the initial conditions given by (4) are chosen as follows: where S 0 = 220870336, E 0 = 16000, I 0 = 4, V 0 = 0 and R 0 = 0. Figures 1 and 2 illustrate profiles of initial conditions. Figure 1 describes initial conditions (3), it has uniform population distribution throughout the domain [−2, 2] , for each sub-populations considered in the study (Fig. 4). It is assumed that the proportion of susceptible individuals is greater than all other sub-populations, while the vaccinated and recovered populations are kept at zero. On the other hand, Fig. 4 illustrates initial condition (4) with susceptible, exposed and infected are concentrated at the origin in the center of the domain [−2, 2] and decrease exponentially on both sides of the origin while the vaccinated and recovered concentrations are considered zeros in this case. (3)

Qualitative analysis
The present section deals with qualitative properties of the model (1), including existence and uniqueness, positivity, and boundedness of the solution to assure its global time existence. Moreover, the stability of the model (1) steady states is discussed in detail.
Well-posedness of model. In the analysis of partial differential equations (PDE), it is important to verify whether the PDE combined with some initial and boundary conditions is well-posed since it claims the solution www.nature.com/scientificreports/ existence and uniqueness. To establish the well-posedness of a mathematical problem under consideration, the semi-group theory approach is employed. In this context, a Banach space denoted by X = C(�; R) is considered. This space consists of real-valued continuous functions φ defined over the set , where X + representing the positive cone of this space. Moreover, the norm |φ| X is defined as the supremum of |φ(x)| over all x ∈ . We assume the existence of a linear operator A ξ , defined as follows: where ξ stands for the respective model state variables, is the differential operator ∂ 2 ∂x 2 and Based on the well-established knowledge pertaining to the operator A ξ , it is known that A ξ serves as the infinitesimal generator of a strongly continuous semigroup denoted by {e tA ξ : t ≥ 0} consisting of linear operators in the Banach space X 14 . Moreover, considering the operator defined as It can be inferred that the above operator also serves as an infinitesimal generator of a strongly continuous semigroup {e tA : t ≥ 0} consisting of linear operators in the Banach space Y = X 514 . Here, the Banach space Y is coupled with a norm defined by: and Y + = X 5 + ⊂ Y indicates a positive cone in the space Y. Further, considering F as a nonlinear operator given Y as: As a result, the model (1) can be rewrite in a comprehensive form as follows: where It is obvious that the operator F is Lipschitz continuous on Y + due the preserving following property where, L = µ + µ 0 + 2(γ + ψ ν + η ν + θ + ω) . Therefore, with the help of Theorem 3.3.3 in 15 , we state the following preposition. (5), and u 0 = (S 0 , E 0 , I 0 , V 0 , R 0 ) T ∈ Y + be the respective initial condition. It is assumed that there exists a unique continuously differentiable solution u(t) of (6) defined on some maximum interval of existence [0, τ ] , satisfying the following equation:

Proposition 3.1 Let A be an operator defined by
Positivity and boundedness. The subsequent lemma, which has been widely employed in demonstrating the non-negativity of the solution to system (1) [16][17][18] .  (7), the expression k(x, t) ∈ C(� × [0, τ )) is a continuous function.
Given Lemma 3.1, we can deduce the following result.
. Then, from the first equation of model (1), we have: It is evident that (8) satisfies the criteria presented in Lemma 3.1. Thus, by applying Lemma 3.1 directly, we obtain S(x, t) ≥ 0 on � × [0, τ ) . Similar arguments can be applied to prove the non-negativity of the remaining state variables E(x, t), I(x, t), V(x, t), and R(x, t).
Next, we demonstrate the boundedness of the solution to model (1) to provide the existence of a global solution. This is accomplished through the subsequent result.

Proof
The Jacobian of the system (1) at the DFE point ξ 0 = ( � d , 0, 0, 0, 0, 0) is, and the associated characteristic polynomial is, where The coefficients of the characteristics equation (9) a 0 , a 1 , a 2 are positive, claimed by R 0 < 1 . It fulfill the necessary condition of Routh-Hurwitz stability criteria 19 . Thus the DFE point is locally asymptotically stable in whenever R 0 < 1 .

Linearized stability of endemic equilibrium.
To discuss linearized stability of endemic equilibrium i.e.
Scientific Reports | (2023) 13:12108 | https://doi.org/10.1038/s41598-023-38925-w www.nature.com/scientificreports/ Global stability. This section presents the global asymptotical stability of the reaction-diffusion system (1) equilibria. The Lyapunov stability theory is used and constructs a suitable Lyapunov function. This theory was presented by a Russian mathematician and until now used as a strong and fundamental tool for investigating the stability of nonlinear dynamical systems 20 . Moreover, this approach does not requires the exact solution of the system examined for stability. We proceed with the following result. (1) is stable globally asymptotically.
Proof In order to establish the result, we define a Lyapunov functional as follows: where and τ indicates the maximumal time.
The time derivative of V(t) is Using condition stated in (2) we have Since, Therefore, for all t ≥ 0 and x ∈ . It follows from system (1) By applying the widely used LaSalle's invariance principle presented in 21 , concluded the disease-free equilibrium ξ 0 of spatio-temporal model (1) is globally asymptotically stable if R 0 < 1 .

Numerical methods
The present section deals with the numerical solution of spatio-temporal epidemic vaccine model (1). For this purpose, two numerical schemes are proposed namely the finite difference operator splitting (FDOSM) and mesh-free scheme based on multi-quadratic radial basis function is take into account. In the first method, the finite-difference approximation is used along with the operator-splitting technique and in the second method radial basis function approximation is implemented with operator-splitting technique. The stability and positivity of these numerical methods for an epidemic spatio-temporal model can be found in 4,6,22 . The aim of the study is to analyze the efficiency and effectiveness of the proposed schemes for epidemic models. The Operator-splitting technique is very useful for solving non-linear spatio-temporal partial differential equations since it effectively Scientific Reports | (2023) 13:12108 | https://doi.org/10.1038/s41598-023-38925-w www.nature.com/scientificreports/ handles the complexity and non-linearity of reaction-diffusion equations. By implementing this technique, the proposed model (1) can be split into two sub-systems, the nonlinear reaction system of equations given by and linear diffusion system of equations, Operator splitting based time descritization. The time descretization of the proposed model (1) can be done in two steps. Initially, using first-order time difference by taking half time step 1 2 dt from 0 to 1 2 dt . The system (16) can be written in descretized form as follow: Next, half time step is taken from 1 2 dt to dt and using first-order time difference, the system (17) is represented in descretized form as: Space descretization. The spatial derivatives in (19) will be approximated in two ways, with the conventional second-order finite-difference and with multi-quadratic radial basis function approximation, describes as follow: Space descretization with finite-difference. The finite-difference approximation for second order derivative is defined as: Space descretization with MQ RBF. In this scenario a set of K centers, Then, the RBF interpolation of a given function f (x) is expressed as: The function φ(r) is termed as radial basis function (RBF), defined for r ≥ 0 where, r = (x i − x j ) 2 + c 2 ≥ 0 , and i, j = 1, 2, 3, . . . , K , while the parameter c is known as shape parameter and not necessary to involve RBF interpolation approximation. Moreover, the coefficients α j in the above expansion can be computed from the following interpolation condition, at a set of nodal points x i for i = 1, 2, 3, . . . , M . In our case, the centers and collocation points are considered similar, utilizing the RBF interpolation yields a given partial differential equation to the following system of linear equations, where α is a K × 1 vector having entries α j for j = 1, 2, . . . , K and B is K × K interpolation matrix with entries, and Moreover, the derivative of a function f (x) can be approximated using RBF as fallow: where and L represent the derivative operator of first or second-order, defined as follow: Scientific Reports | (2023) 13:12108 | https://doi.org/10.1038/s41598-023-38925-w www.nature.com/scientificreports/ The symbols represent the interior and ∂� denotes the boundary of the domain of study. Now, the proposed procedure based on MQ RBF known as a meshless method, is applied to solve the spatio-temporal vaccine model (1). It has been proven that this method have significant approximation capabilities. The attractive feature of the meshless method is that it is easily extendable to higher dimensions in both cases on scattered and uniform data. By utilizing the RBF approximation for the spatial derivatives in (19), the system of linear equations (19) for the half time step 1 2 dt to dt can be expressed as: where and F is a matrix with entries f i , defined below: where ξ indicates S, E, I, V and R and, The entries of matrix B d are The coefficients α n+ 1 2 rj for r = 1, 2, 3, 4, 5 and j = 1, 2, 3, . . . , K are computed from (21). The system of equations (21) can be solved either by using standard inversion scheme, LU-factorization or by Gauss-elimination technique.

Simulation and discussion
The reaction-diffusion COVID-19 vaccine model (1) is simulated by using the proposed iterative schemes FDOSM and RBF as discussed previously. The simulation results are performed in Matlab version R2022b for the time period 0-700 days. The coefficients of diffusivity are assumed as α S = 0.00005, α E = 0.0005, α I = 0.001, α V = 0.001 and α R = 0 . Initially, the model (1) is simulated for the values of the embedded parameters provided in Table 1 with and without diffusion, in order to analyze the role of diffusion on the dynamic of the COVID-19 epidemic. In the case of without diffusion, the coefficients of diffusivity are considered as α S = 0, α E = 0, α I = 0, α V = 0, α R = 0 . Figure 3 describes the dynamics of symptomatic and exposed individuals with and without diffusion. A notable observation is that the presence of diffusion leads to a significant reduction in the number of infected individuals compared to the scenario without diffusion. According to biological facts, diffusion phenomena help to restrict the public gathering, as a result, it yields lower chances of getting the infection.

Simulation using FDOSM.
In this subsection, model (1) simulation is performed with and without diffusion by utilizing a finite difference operator-splitting approach. For this purpose, the time step is taken t = 0.028 and corresponds to the number of spatial points N = 60 . The spatial step is x = 0.06 which is chosen based on Von Neumann stability criteria presented in 23 . The simulation is carried out for both uniform and nonuniform initial conditions given in Eqs. (3) and (4) respectively.
Simulation based on uniform initial condition (3). This part of the paper presents the dynamics of exposed and symptomatically infected populations based on uniform initial conditions (3). The simulation is performed with and without diffusion at x = 0 and x = 1 . The resulting plot is depicted in Fig. (4) where the subplots (a) and (b) demonstrate the dynamical aspects of exposed and infected classes at x = 0.0 , and the subplots (c) and (d) present the population dynamics at x = 1.0 . The values of parameters used in the simulation are given in Table 1. The population is a spatially uniformly distributed due to the initial condition (3). Therefore, one can observe similar behavior with and without diffusion at x = 0.0 and as well as at x = 1.0.
Simulation based on non-uniform initial condition at x = 0. This subsection presents the simulation results describing the time evolutionary trajectories of the model (1) with and without diffusion at x = 0 using FDOSM. Moreover, this section illustrates the role of various parameters including β , κ , ψ ν and η ν with diffusive and non-diffusive cases. These parameters correspond to the implementation of various pharmaceutical and nonpharmaceutical interventions such as social distancing measures and vaccination etc. These results are dem-   Figure 5 illustrates the impact of control measure β (the effective contact rate) on the exposed and infectious population. The dynamics are studied in both diffusive and non-diffusive cases as shown in Fig. 5a-d. The profiles of the exposed and infected populations are obtained for the baseline value β = 0.4710 given in Table 1 and it is decreased by 10%, 20% and 30% . It is observed that a reasonable reduction is analyzed in the aforementioned population classes. Moreover, the lowest projected peak up to 97% reduction in exposed and symptomatically infected individuals is noticed with 30% decrease in social contacts in presence of diffusion. The detailed description is presented in Table 2. Further, we demonstrate the impact of parameter κ (the relative transmissibility rate of infection due to exposed individuals) population dynamics in both diffusive and non-diffusive cases as shown in Fig. 6a-d. The graphical results are analyzed initially for the baseline value, i.e. κ = 0.3378 in case of diffusion and without diffusion. One can observe the difference between the projected peak value presented in Table 2. Furthermore, the value of κ is reduced by 10%, 20% , and 30% respectively and as a result, 92% decrease is observed with diffusion. Figure 6 describes the role of the implementation of vaccination control measure ψ ν on the disease incidence. The simulation is performed for the baseline value of this parameter i.e., ψ ν = 0.0209 . Then ψ ν is enhanced to 10%, 20%, 30% and 40% to the baseline value. It is noticed that exposed and symptomatically infected individuals get decreases and with 40% enhancement in vaccination rate, a 95% decrease is observed in the case of diffusion. According to biological facts, vaccines contain anti-viral agents that are responsible for boots the immunity of suspects and it will probably increase the recovery rate. That's why a significant reduction in the infected population is observed. Figure 8 depicts the role of vaccine waning rate η ν on exposed and symptomatically infected individuals. The behavior is analyzed by decreasing 10%, 20%, 30% and 40% vaccine waning rate from the baseline value η ν = 0.0132 . A significant reduction is analyzed with 40% decrease in vaccine waning rate. The more effective result is observed in the case of diffusion.
Simulation based on nonuniform initial condition at x = 1.0. The current section presents the dynamics of exposed and infected individuals based on nonuniform initial conditions (4) at x = 1.0 using FDOSM. The simulation is performed for the model in diffusive and non-diffusive cases. The dynamics of exposed and infected populations with and without diffusion under the variation in β is examined in Fig. 9. According to the initial conditions profiles given by Fig. 4, the population concentration is low at x = 1.0 . Therefore, in the absence of diffusion, the exposed and infected individuals are decreasing exponentially by decreasing the control measure β . On the other hand, when diffusion is involved the exposed and infected individuals rise for the baseline value β = 0.4710 . Moreover, with the reduction in β by 10%, 20% and 30% to the baseline value, a reasonable reduction is observed in exposed and infected individuals. The projected peak values, in this case, are presented in Table 3. Figure 10 depicts the time evolutionary trajectories obtained by using FDOSM for the control measure κ with and without diffusion at x = 1.0 . According to Fig. 4, the exposed and infected individuals have low concentrations at x = 1.0 . Therefore, in the absence of diffusion, the number of exposed and infected individuals decreased exponentially and no significant reduction is observed by reducing κ . Moreover, in the case of diffusion i.e. the population moves from higher concentration to lower concentration, the exposed and infected population increased initially for the baseline value at x = 1.0 . Further, the simulation is performed for 10%, 20% and 30% reduction in κ where the lowest peak is noticed for exposed and infected individuals with 30% reduction in the value of κ as shown in Table 3. Exposed Individuals E(t) 10 7 Dynamics without diffusion Dynamics with diffusion (a) Infected Individuals I(t) 10 6 Dynamics without diffusion Dynamics with diffusion (b) Figure 3. Impact of diffusion on the exposed and infectious individuals (a) shows profile of the exposed population and subplot (b) depicts the dynamics of infectious individuals with and without diffusion.
Scientific Reports | (2023) 13:12108 | https://doi.org/10.1038/s41598-023-38925-w www.nature.com/scientificreports/ One of the most important control interventions against the COVID-19 pandemic is proper vaccination. The impact parameter describing the vaccination rate ψ ν on the exposed and infected population dynamics at spatial point x = 1.0 with and without diffusion is shown in Fig. 11. The evolutionary curve is obtained initially for the baseline value of ψ ν i.e., 0.0209. Furthermore, the vaccination rate is enhanced by 10%, 20%, 30% and 40% to the actual value in order to analyze the dynamics of the infected population. According to initial condition (4) the population concentration at x = 1.0 is low. As a result, if diffusion is not considered i.e., there is no spatial movement of the population, the exposed and infected population is decreased exponentially and enhancing the vaccination rate has no reasonable impact on the reduction of infection. On another hand, in the case of diffusion, as the population moves from higher concentration to lower concentration, the number of infected individuals initially increases for the baseline value of ψ ν . By enhancing the vaccination one can observe a significant reduction in the infected population. Moreover, the lowest peak is observed in the case of 40% reduction in ψ ν with diffusion.
The impact of vaccine waning rate η ν on the respective population dynamics is described in Fig. 12 with and without diffusion at x = 1.0 . The behavior is observed for the baseline value 0.0132 of η ν and then it is reduced by 10%, 20%, 30%, 40% , 50% . From the evolutionary curve obtained by the FDOSM scheme, it is observed that in the absence of diffusion, the number of infected and exposed individuals reduces exponentially due to the low concentration of population at x = 1.0 . In the case of diffusion, the number of infected individuals increases   Table 4. It is noticed that the results obtained from both numerical techniques are similar. But the benefit of using mesh-free operator-splitting based on multi-quadratic radial basis function is that this approach is a unique distinction as it works on scattered data and does not require underlying meshes or structured nodes. This approach is easily extendable to higher dimensional reaction-diffusion models. Moreover, it produces a better result on a low number of collocation points as compared to the finite-difference operator splitting scheme.
Simulation based on initial condition (4) at x = 1.0. Finally, in Figs. 21 and 22, we present mesh plots of model (1) to study the combined spatial and temporal dynamics of different populations on the entire time and spatial points in [−2, 2] . The mesh plots associated with the proposed operator-splitting meshless numerical scheme demonstrate a consistent behavior. It is evident that the scheme effectively maintains the positivity property of the solution, ensuring that the solution remains positive for all t > 0 and spatial points within the defined domain [−2, 2] . Moreover, the solution converges toward the steady states. In addition, according to the initial population distribution given by (4) concentration of the population is higher at the center of the spatial domain for the time period 0-10 days and then gradually diffuses in [−2, 2] with time. Therefore, the exposed population is increasing due to the spatial movement of susceptible and as a result of interaction with infected and exposed individuals. Therefore, the population of infected as well as the vaccinated compartments increased.

Conclusion
In the present study, a spatio-temporal COVID-19 vaccine model is formulated and analyzed qualitatively and numerically. The aim of this work is to analyze the impact of diffusion as well as some control interventions including social distancing (by reducing effective contacts), vaccination rate and vaccine waning rate on the dynamics of the infected population in a spatially heterogeneous environment. Initially, the spatio-temporal COVID-19 model is analyzed qualitatively in detail. The basic mathematical analysis including the existence, uniqueness, boundedness and positivity of the solution is presented. The local and global stability of the spatio-temporal model at the steady sates is shown using well-known techniques. Moreover, the model is solved Exposed Individuals E(t) 10 6 Dynamics with diffuison  www.nature.com/scientificreports/ numerically based on uniform and non-uniform initial conditions with two different numerical schemes named: finite difference operator-splitting and mesh-free operator-splitting based on multi-quadratic radial basis functions. A detailed simulation is presented using both numerical schemes for diffusive and non-diffusive cases at two different spatial points x = 0 and x = 1 . Two different spatial points are chosen due to the fact that x = 0 represents the highly populated area as compared to x = 1 . The simulation results concluded that based on the initial conditions (3) and (4) at different spatial locations, the diffusion effect will play a significant role in curtailing infection transmission in a heterogeneous environment with concentrated regions. Further, it is observed that reducing β to 30% in the case of diffusion, the infection will be reduced to 97% while enhancing vaccination rate to 40% a 95% decrease is observed in the infection incidence. The restriction of public gatherings during COVID-19 outbreak is important as the virus is airborne and easily transmit with low social distancing. Thus, it is concluded that implementing the suggested control strategies with diffusion is beneficial and helps in the eradication of infection from the community. Exposed Individuals E(t) 10 6 Dynamics with diffuison (d) Figure 6. Impact of effective contact rate κ , with and without diffusion at x = 0.0 on exposed individuals shown in subplots (a,b) and on infected individuals shown in subplots (c,d).  Exposed Individuals E(t) 10 6 Dynamics with diffuison Exposed Individuals E(t) 10 6 Dynamics with diffuison Exposed Individuals E(t) 10 4 Dynamics with diffuison   Exposed Individuals E(t) 10 4 Dynamics with diffusion Exposed Individuals E(t) 10 4 Dynamics with diffusion Exposed Individuals E(t) 10 4 Dynamics with diffuison  Exposed Individuals E(t) 10 6 Dynamics with diffuison (d) Figure 13. Impact of effective contact rate β , at x = 0.0 using RBF. The subplots (a,b) illustrate the dynamics of exposed individuals with and without diffusion while the subplots (c,d) present the respective profile of infected individuals with and without diffusion. Exposed Individuals E(t) 10 6 Dynamics with diffuison (d) Figure 14. Impact of transmissibility rate κ at x = 0.0 using RBF. The subplots (a,b) illustrate the dynamics of exposed individuals with and without diffusion, and the subplots (c,d) demonstrate the respective profiles of infected individuals with and without diffusion. Exposed Individuals E(t) 10 6 Dynamics with diffuison (d) Figure 15. Impact of vaccination rate ψ ν , at x = 0 using RBF. The subplots (a,b) illustrate the dynamics of exposed individuals with and without diffusion while the subplots (c,d) show the respective profiles of infected individuals with and without diffusion. Exposed Individuals E(t) 10 6 Dynamics with diffuison Exposed Individuals E(t) 10 4 Dynamics with diffusion (d) Figure 17. Impact of effective contact rate β at x = 1.0 using RBF where the subplots (a,b) illustrate the dynamics of exposed individuals with and without diffusion and the subplots (c,d) show the respective profile of infected individuals with and without diffusion.  Exposed Individuals E(t) 10 4 Dynamics with diffusion (d) Figure 18. Impact of transmissibility rate κ at x = 1.0 using RBF. The subplots (a,b) illustrate the dynamics of exposed individuals with and without diffusion and the subplots (c,d) analyze the respective profiles of infected individuals with and without diffusion. Exposed Individuals E(t) 10 2 Dynamics without diffusion Exposed Individuals E(t) 10 4 Dynamics with diffusion Exposed Individuals E(t) 10 4 Dynamics with diffusion